Magnetic and nematic orderings in spin-1 antiferromagnets with single-ion anisotropy 
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We study a spin-one Heisenberg model with exchange interaction, J, uniaxial single-ion exchange 
anisotropy, D, and Zeeman coupling to a magnetic field, B, parallel to the symmetry axis. We com- 
pute the (D I J, B I J) quantum phase diagram for square and simple cubic lattices by combining an- 
alytical and Quantum Monte Carlo approaches, and find a transition between XY-antiferromagnetic 
and ferronematic phases that spontaneously break the U(l) symmetry of the model. In the language 
of bosonic gases, this is a transition between a Bose-Einstein condensate (BEC) of single bosons and 
a BEC of pairs. Our work opens up new avenues for measuring this transition in real magnets. 
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The unambiguous realization of BECs in laser cooled 
collections of cold atoms [I] [2] triggered the search for 
more exotic states of matter and phase transitions that 
take place in bosonic gases. In parallel with these efforts, 
other experimental groups demonstrated that several as- 
pects of the BEC quantum phase transition can also be 
measured in quantum magnets that are alternative re- 
alizations of bosonic gases [3H7]. The main advantage 
of this alternative approach is that the much lighter bo- 
son mass leads to much higher transition temperatures, 
while the uniform character and well-defined temperature 
of quantum magnets are crucial for the study of quan- 
tum phase transitions. The main disadvantage is that 
the continuous symmetry associated with particle num- 
ber conservation in atomic gases is always an idealization 
for the case of quantum magnets [8]. However, many real 
magnets, for which the continuous symmetry breaking 
terms are much smaller than the ordering temperature, 
allow for measuring the universal behavior of continuous 
symmetry breaking critical points over a large window 
of temperatures. The simplest example is the magnetic 
field induced BEC quantum critical point (QCP) that is 
observed in several quantum magnets |3HZ]. 

5=1 magnets can be mapped into gases of semi- 
hard core bosons via a generalization of the Matsubara- 
Matsuda transformation [9j [10] that also maps the local 
magnetization into the boson density: nj = 5| + 1. In 
contrast to the hard-core bosons associated with S = 1/2 
magnets, it is possible to study "Hubbard-like" bosonic 
gases with on-site density-density interactions because 
nj < 2 [TTHT5] . Moreover, the semi hard-core constraint, 
rij < 2, can be incorporated as an infinitely repulsive 
three-body on-site term that precludes phase segrega- 
tion when the two-body term is attractive. This situa- 
tion is ideal for studying transitions between single-boson 
BECs and BECs of pairs, whose counterparts in atomic 
physics are transitions between BECs of atoms and di- 
atomic molecules (transitions between XY-magnetic and 
nematic orderings in the spin language). 

While the XY terms of the exchange interaction play 



the role of the kinetic energy, the Ising terms map into 
off-site density-density interactions. On-site density- 
density interactions are generated by uniaxial single-ion 
anisotropy terms of the form D(Sj) 2 . A magnetic field, 
B, parallel to the symmetry axis acts as a chemical 
potential because it couples to the total magnetization 
M z = J2j Sj that coincides with the total number of 
bosons up to an irrelevant constant. Previous works have 
exploited this spin-boson mapping for studying spin su- 
persolid phases of S = 1 Heisenberg models with uniaxial 
exchange and single- ion anisotropics [T6HT8] . The main 
purpose of this work is to study the quantum phase di- 
agram (QPD) of the isotropic S = 1 Heisenberg model 
when the single-ion anisotropy and Zeeman terms are in- 
cluded. This model is relevant for describing several Ni- 
based organic compounds as well as inorganic systems 
that are discussed at the end of this work. 

Our 5 = 1 model is defined on a hypercubic lattice, 

H = jJ2Si-S j + Dj2(S!) 2 -Bj2st, (1) 

and the antiferromagnetic (AFM) exchange coupling J > 
only connects nearest-neighbor sites Since an 

attractive on-site interaction is needed for pairing the 
bosons, we will only consider the D < case that corre- 
sponds to easy-axis anisotropy. 

Hamer et at. computed the mean field QPD of H on a 
square lattice and obtained a single phase transition from 
AFM Neel (AFM-z) to the fully polarized (FP) state for 
large \D\/J values [19]. However, an intermediate ne- 
matic phase must exist in this regime according to the ef- 
fective low-energy model, that is derived by expanding 
in the small parameter J/\D\ (strong coupling expansion 
in the bosonic language). The low-energy subspace, 5, is 
a direct product of the S z = ±1 doublets of each lattice 
site that are separated from the S z = states by an en- 
ergy gap \D\ [20 . H is obtained by applying a canonical 
transformation and projecting into S: H = Pse K T~Le~ K Ps 
(k is the anti-hermitian generator of the canonical trans- 
formation and Ps is the projector on S). If we use a 



2 



pseudospin 5 = 1/2 variable to describe each doublet, 
S z = 2s z , we obtain the following expression for ti up to 
quadratic order in J: 



H 



(2) 



with J z = AJ-J 2 /D and = j^/jj ft is an s = l / 2 
XXZ model whose QPD is well-known in any dimension. 
In particular, the mean field phase diagram is qualita- 
tively and quantitatively correct for spatial dimension 
d > 2. Since the case of interest corresponds to effective 
strong easy- axis anisotropy, \J Z / J xy \ ~ 4\D\/J ^> 1, 
the low field ground state has Neel ordering that extends 
up to the spin-flop field, B s ^ whose mean field value is 
B s f ~ dy/(J z + J x y)(J z - J*y)/2. The corresponding 
curve is the lower dotted line on the left of Figs. [IJi and 
b. At B = £? s f, the pseudospin variables flop from the 
Neel state to a canted ferromagnetic state (J xy < 0) 
whose canting angle relative to the z-axis is given by 
cosa s f = 2(sj) 



^(J z + J x y)/(J Z -J x y). The effec- 



tive operator for Sj~ is Sj~ 



P S e K S+ e 



This identity implies that the planar ferromagnetic order- 
ing in the pseudospin variables or "spin-flop" phase cor- 
responds to ferronematic (FNM) ordering in the original 
S = 1 spin variables, i.e., (st) ^ implies (Sj~ ) ^ 0. On 
the other hand, the effective operator for is equal to 



zero, Sj~ 



Pse K Sj~e K Ps = 0, because of the following 



symmetry argument. Being odd under time reversal sym- 
metry, Sj~ must be equal to a polynomial form that only 
contains odd terms in the s\ variables v = {x, y, z). Such 
polynomial form must be odd under a 7r spin rotation 
along the z-axis because e t7V ^ s iSj~e~ l7r ^ s * = —Sj~. 
Since e l7r ^ s is~j~e~ l7r ^ s * = si~, the polynomial form 
must be equal to zero, implying absence of planar mag- 
netic ordering in the large D/J limit and confirming the 
FNM character of the intermediate phase. The transi- 
tion to the fully saturated state is of second order in this 
regime and belongs to the BEC universality class in di- 
mension d+2. A mean field treatment [19 of the original 
Hamiltonian, %, misses the second order fluctuations in 
J (J xy = J 2 jD) that stabilize the FNM phase. 
The approximated value of the saturation field is 



B sat (\D\/J»l)^^(J z 



JX 



d(2J 



11 
D 



), (3) 



and the corresponding curve is the upper dotted line 
on the left of Figs. [IJi and b. While Eq. ([3| is a good 
approximation for B sat if \D\/J ^> 1, the exact curve 
Bsat(\D\/J ^> 1) can be computed analytically as long 
as the transition remains continuous. By solving the two 
body problem of diagonalizing T~L in the M z = N — 2 
invariant subspace (two flipped spins relative to the FP 
state) we obtain the exact energy, E g (M z = N — 2), of 
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FIG. 1: Quantum phase diagrams in (a) square lattice and 
(b) cubic lattice. Double and single solid lines correspond to 
first and second order phase transitions respectively. Single 
solid lines are obtained by exact solutions while double solid 
lines are guides for the eye based on QMC results. Insets: 
enlarged views of the FNM to AFM-xy transition. Dashed 
and dotted lines are analytical solutions for the weak and 
strong anisotropy regimes respectively. 



the two bosons bound sate. The condensation of these 
pairs leads to the FNM ordered state. If the transition is 
continuous, the exact value of B sat is the field such that 
E g (M z = N-2) = E g (M z = N)= N(Jd + D-B). The 
resulting curve is shown as a full line in the \D\/J ^> 1 
region of Figs. |TJl and b. 

The opposite limit, \D\/J <C 1, is well described by 
a mean field treatment of the original Hamiltonian T~L. 
The mean field Neel state |^n) = l-Oj I-Oj * s 

the most stable at low-enough fields. Here A and B de- 
note the two sublattices, while |l)j, |l)j and |0)j are the 
eigenvectors of £J with eigenvalues 1, —1 and respec- 
tively. In this regime there is a spin-flop transition, but 
to a canted XY AFM (AFM-xy) phase that is described 
by the mean field state 



I* 



sf/ 



e iQ ' r i sin0[cos0|l)j + sin^l)^ +cos<9|0) j , 



where Q is the AFM wave vector that has all the com- 
ponents equal to tt. The optimal variational parameters, 
0o and 0o 5 are obtained by minimizing the mean field en- 
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ergy (ty B { |%|\I/ s f). The dashed line on the right of Figs.[TJi 
and b corresponds to the spin flop curve B S {(D/J) that 

results from <* N |ft|* N ) = <* B f(0o,0o)|«|*Bf(0o,&>)>. 

For small D/J, the second order transition between 
the spin-flop and FP states belongs to the BEC univer- 
sality class. The exact value of the saturation field in this 
regime is shown as a full line in the upper right region of 
Figs, [l^i and b, and given by the equation 



B = D + 4dJ. 



(4) 



However, we know that the effective interaction between 
bosons should become attractive for \D\ > \D C \\. There- 
fore, the second order transition line described by Eq.Q 
should become of first order at a tricritical point (TCP) 
with coordinates [D c \/ J,Bq(D c \/ J) / J\ (see Figs.JT^i and 
b). The region near the TCP is well described by the 
Ginzburg-Landau (GL) free energy density, 



f(<j>) = (B-B )\<j ) \ 2 +u\<j ) \ i +wl 



(5) 



Here <ft is the complex order parameter for the BEC of 
single bosons, u and w are the amplitudes of the ef- 
fective two-body and three-body interactions in the long 
wavelength (or continuum) limit. The field induced tran- 
sition is continuous for repulsive u > and it happens at 
B sat = Bq [see Eq.Q]. However, it is clear from Eq. ([5| 
that the transition becomes discontinuous for u < 0. 
In this case, the transition field is B sat = Bq + u 2 /Aw 
and the discontinuous change of the boson density is 
Am z = A |^| 2 = -u/2w (m z = T,j(Sj)/ N )- Th e am- 
plitude u changes sign when the two-boson scattering 
length, a s , diverges in d = 2 (a s — » oo) and becomes 
equal to zero in d = 3 (a s = 0). These conditions deter- 
mine the values of D c i / J in d = 2 and d = 3 respectively, 
that can be obtained by computing the effective interac- 
tion vertex in the long wave length and low frequency 
limits: 



T q (k,k';u) = V q (k)- 



d d p V q - p (k)T p (k,k f ;uj) 



n (2n) d uj - e k+p - e k > 



-iS 
(6) 



where V q (k) = 2D + j q + (y/2 - 2)( 7fc+q + 7fc ), e q = 
(2d J + 7 q ) is the single boson dispersion at the TCP 
and 7 q = 2J^2 v cosk u . By solving Eq. ( |6| ) for q — 0, 
k = k' = Q. and uj — » 0, we obtain \D C \\/J = 4d/3. 
A similar analysis cannot be applied to the point where 
the FP phase boundary changes from second to first or- 
der coming from the strongly anisotropic side \D\ ^> J. 
Note that the FNM phase disappears right at this criti- 
cal D = D C 2/J point (see Fig. [I]), while the magnetiza- 
tion vs field curve becomes discontinuous (see insets of 
Figs. [2] and [3|. This discontinuity indicates that it is a 
critical end point (CEP) and the effective GL theory is 
not applicable. Then, the coordinates of the CEP must 
be obtained from the quantum Monte Carlo (QMC) sim- 
ulations that we describe below. 
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FIG. 2: Magnetic field sweep showing AFM-z to FNM to 
FP sequence of ground states for d — 2 lattices of L 2 sites 
and periodic boundary conditions (PBC). Data shown are for 
D/J — —5 and inverse temperature of /3J — 4L. The in- 
set shows m z along the line B(D/J) where the two-magnon 
ground state is degenerate with the FP state. 



Our analysis of the two opposite regimes, \D\/J ^> 1 
and \D\/J < 1, indicates that there is a transition be- 
tween AFM-xy and FNM orderings in the intermediate 
region. We use a QMC method with global updates [21] 
for studying this regime, because there is no small pa- 
rameter for validating an analytical approach. Although 
H does not have a negative sign problem, standard QMC 
algorithms cannot output the off-diagonal FNM correla- 
tor because of a slowing down problem. Therefore, we use 
a novel multi-discontinuity algorithm [22], that is based 
on the directed-loop algorithm [23] and eliminates the 
problem. 

The different phases are characterized by computing 
the zero frequency AFM-xy and FNM susceptibilities, 

fk ^ j\st(T)S+(T)ST(p)ST(p))dT, 



FNM 



where N = L d is the number of lattice sites. We also 
compute standard thermodynamic quantities, like the 
magnetization, m 2 , and the spin stiffness, p s (response 
of the system to a twist in the boundary conditions), 
that is obtained from the fluctuations of the world lines 
winding numbers along the principal axes [24] . 

Figs. [T^i and b include the d — 2 and d = 3 QPDs 
obtained from our QMC results. Figs. [2] shows the four 
different observables computed as a function of B/J for 
D/J = —5 and different system sizes. Except for the 
FNM-AFM-xy transition, the phase boundaries of the 
first order phase transitions, shown in Figs.[TJi and b, are 
determined from the size dependence of the discontinu- 
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FIG. 3: QMC results for d = 3 lattices of L 3 spins with 
PBC. The temperature is scaled with the systems size, f3J — 
4L, and the Hamiltonian parameters are (a) D/J — —5, (b) 
D/J = —7.5. Inset of (b) shows m z along the line B(D/J) 
where the two-magnon ground state is degenerate with the 
FP state (PJ= 6L). 




FIG. 4: Finite size scaling of the QMC results for the AFM- 
xy structure factor, S + ~ (Q) = S+ e iQ<r i~ ri) /N , near 
the FNM to AFM-xy transition in d = 2 with B/J — 4.1 and 
f3 — 4L. The values of the exponents for the Ising universality 
class in dimension 3 = 2 + 1 are extracted from Ref. 25 . The 
common crossing point at D c / J — —4.22 ± 0.02 (see inset), 
together with the curve collapse in the critical region, suggests 
that the transition is continuous. 



ity of the uniform magnetization and the corresponding 
kink in the energy density. These boundaries agree very 
well with our analytical solutions in the limiting regimes 
\D\/J > 1 and \D\/J < 1. Moreover, the QMC re- 
sults indicate that the transition to the saturated state 
becomes of first order between the CEP and TCP that 
we discussed above. The coordinates of the TCP coin- 
cide with the exact values obtained by solving Eq.|6]). 
The coordinates of the CEP obtained from our QMC 
simulations are D C 2/J — —3.70 ± 0.03 for d = 2 and 
D c2 /J = -6.33 ± 0.03 for d = 3. Figs. [§i and b show 
that the magnetization curve has a small discontinuity 
for D/J = —5, while it is continuous for D/J = —7.5. 
The first order transition line for |D c i| < \D\ < \D c2 \ 
falls consistently above the curves given by Eqs. ([3| and 
Q, as expected from our GL analysis near the TCP. 

The transition from FNM to AFM-xy ordering (double 
full-dashed line in Figs, [l^i and b) spontaneously breaks 
the discrete symmetry of global spin rotation by tt along 
the z-axis. Consequently, if continuous, this transition 
should belong to the Ising universality class in dimension 
d+ 1. The scaling analysis shown in Fig. [4] indicates that 
this transition is most likely continuous away from the 
FP and AFM-z phases. However, our magnetization vs. 
field curves indicate that it becomes weakly first order 
near the boundaries with these two phases, implying that 
the upper end of the FNM to AFM-xy phase boundary is 
a CEP, while the lower end corresponds to a triple point 
(TP) at the junction of the FNM, AFM-xy, and AFM-z 
phases (see Fig. [T]). 

The continuous or quasi-continuous nature of the FNM 
to AFM-xy quantum phase transition indicates that the 
single-boson condensate is continuously converted into 



a condensate of pairs (the condensate density is equal 
to the particle density p = 1 — m z in the low-density 
limit). This observation implies that BECs of pairs and 
single-bosons coexist in a finite region of the AFM-xy 
phase that ends up at the phase boundary between the 
two phases where the single-boson BEC disappears com- 
pletely: (S+) = 0. Indeed, for d = 3 and D = -6.3J, 
the size of the boson-pair, £ 0.77, is three times shorter 
than the average inter-boson distance, p -1 ^ 3 , right below 
the saturation field (p = 1 — m z ~ 0.1). 

The AFM-xy and FNM orderings correspond to BECs 
of single bosons and pairs of bosons, respectively. The 
shape of the phase boundary opens the possibility of mea- 
suring magnetic field induced transitions between these 
two phases (see Figs. [IJi and b). Since a direct experi- 
mental detection of the spin-nematic order parameter can 
be rather challenging, our predictions for the quantum 
phase diagram and behavior of different thermodynamic 
properties are crucial for unveiling this ordering in real 
magnets. While many 5 = 1 magnets are described by T~L 
[26H29] , it is vital to know what are the optimal ratios of 
D/J for detecting the FNM ordering and characterizing 
the FNM to AFM-xy quantum phase transition. Since 
most of these compounds are organic magnets, the D/J 
ratio can be largely tuned as a function of pressure. Thus, 
knowing the appropriate range of D/J values is necessary 
for detecting organic materials in which such a transition 
can be induced by pressure in magnetic fields that should 
be nearly 95% of the saturation field. 

Finally, we mention that field-induced spin supersolid 
states (coexistence of AFM-z and FNM orderings) exist 
at least in the strongly anisotropic limit of T~L for trian- 
gular [30H32] and face-centered-cubic lattices [33]. Other 
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exotic states have been reported for the kagome lattice 
[20] . Ferronematic order has also been obtained for S = 1 
Heisenberg models that include biquadratic interactions 

[TO1I34H37]. 
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